Contexto

• Sistemas em paralelo aumentam a confiabilidade
• Redundância em setores críticos:
  • Energia elétrica
  • Servidores de internet
    
  • Sistemas de transporte
  • Processos industriais

Objetivo do Projeto

• Modelar sistema com 2 componentes idênticos em paralelo
• Analisar falhas (p) e reparos (q) como cadeia de Markov
• Simular dinâmica do sistema em R
• Estudar confiabilidade do sistema

Estados do Sistema

• Estado 2: Ambos componentes funcionando
• Estado 1: Um funcionando, um falhado
  
• Estado 0: Ambos falhados

Probabilidades de Transição

• p: Probabilidade de falha (componente funcionando)
• q: Probabilidade de reparo (componente falhado)

Cálculo das Probabilidades

Do estado 2: - Permanecer em 2: (1-p)² - Ir para 1: 2p(1-p) - Ir para 0: p²

Do estado 1: - Ir para 2: (1-p)q - Permanecer em 1: pq + (1-p)(1-q) - Ir para 0: p(1-q)

Do estado 0: - Ir para 2: q² - Ir para 1: 2q(1-q) - Permanecer em 0: (1-q)²

Matriz Resultante

\[ P = \begin{pmatrix} (1-q)^2 & 2q(1-q) & q²\\ p(1-q) & pq + (1-p)(1-q) & (1-p)q\\ p² & 2p(1-p) & (1-p)² \end{pmatrix} \]

Verificando a validade da matriz

• Uma matriz de transição é válida se satisfaz as duas propriedades fundamentais:

  1. \(P(x,y) \geq 0, \hspace{5mm} \forall x,y \in \mathcal{S}\).

  2. \(\displaystyle \sum_{y \in \mathcal{S}}P(x,y) = 1, \hspace{5mm} x \in \mathcal{S}\).

Verificando a validade da matriz

• Soma na linha 1 da matriz: \[\begin{align*} (1-q)^2 + 2q(1-q) + q^2 &= (1-q)[(1-q) + 2q] + q^2\\ &= (1-q)(1+q) + q^2\\ &= 1 - q^2 + q^2\\ &= 1. \end{align*}\]

Verificando a validade da matriz

• Soma na linha 3 da matriz: \[\begin{align*} p^2 + 2p(1-p) + (1-p)^2 &= p^2 + 2p - 2p^2 + 1 - 2p + p^2\\ &= 2p^2 - 2p^2 + 2p - 2p + 1\\ &= 1. \end{align*}\]

Classificação dos estados

Determinando a distribuição estacionária

• Uma cadeia de Markov de recorrência positiva e irredutível (como a nossa) possui uma única distribuição estacionária \(\pi\).

• Além disso, satisfeitas as condições \(0 < p < 1\) e \(0 < q < 1\) a cadeia é aperiódica,e portanto ergódica. Logo, a distribuição estacionária existe e é igual a distribuição limite da cadeia.

• Podemos obter essa distribuição em função de \(p\) e \(q\) ao resolver o sistema de equações resultantes de \(\pi = \pi P\).

Determinando a distribuição estacionária

\[ \begin{cases} (1 - q)^{2}\pi_{0} + p(1 - q)\pi_{1} + p^{2}\pi_{2} = \pi_{0} \\ 2q(1 - q)\pi_{0} + (pq + (1 - p)(1 - q))\pi_{1} + 2p(1 - p)\pi_{2} = \pi_{1} \\ q^{2}\pi_{0} + (1 - p)\pi_{1} + (1 - p)^{2}\pi_{2} = \pi_{2} \\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]

Determinando a distribuição estacionária

• De L1 e L3, respectivamente:

\[ \begin{cases} \pi_{0}​[1−(1−q)^{2}] = \pi_{1}​p(1−q) + \pi_{2}p^{2}\\ \pi_{2}​[1−(1−p)^{2}] = \pi_{0}​q^{2} + \pi_{1}​(1−p)q \\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]

Determinando a distribuição estacionária

• Arranjando os termos:

\[ \begin{cases} \pi_{0}​[1−(1−q)^{2}] = \pi_{1}​p(1−q) + \pi_{2}p^{2}\\ - \pi_{0}​q^{2} = \pi_{1}​(1−p)q - \pi_{2}​[1−(1−p)^{2}]\\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]

Determinando a distribuição estacionária

• Defina:

\[ \begin{cases} A = 1−(1−q)^{2} = 2q−q^{2} \\ B = p(1−q) \\ C = p^{2} \\ D = 1−(1−p)^{2} = 2p−p^{2} \\ E = q(1-p) \\ F = q^{2} \end{cases} \]

Determinando a distribuição estacionária

• Tem-se então o sistema:

\[ \begin{cases} B\pi_{1} + C\pi_{2} = A\pi_{0} \\ E\pi_{1} - D\pi_{2} = - F\pi_{0} \\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]

Determinando a distribuição estacionária

• Um caminho para resolver o sistema seria multiplicando-se L1 por D, L2 por C e somando as equações. Com o auxilio de L3, por fim, obtém-se:

\[ \begin{cases} \pi_{0} = \frac{p^{2}}{(p+q)^{2}} \\[6pt] \pi_{1} = \frac{2q}{p} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{2pq}{(p+q)^{2}} \\[6pt] \pi_{2} = \frac{q^{2}}{p^{2}} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{q^{2}}{(p+q)^{2}} \end{cases} \]

Determinando a distribuição estacionária

• Observa-se que \(\pi_{0} + \pi_{1} + \pi_{2} = \frac{p² + 2pq + q²}{(p+q)²} = \frac{(p+q)²}{(p+q)²} = 1\)

Simulando a cadeia

• Agora vamos realizar a implementação de uma simulação da cadeia desse sistema no R.

• Partindo do estado 2, simularemos 1000 etapas da cadeia e iremos comparar as proporções dos tempos em que cadeia permance em cada estado com a distribuição estacionária.

Escolha de p e q

Para a simulação que vem logo a seguir, os valores de \(p\) e \(q\) são escolhidos de forma arbitrária considerando os seguintes pontos:

p <- 1 - 0.65
q <- 1 - 0.35

• A probabilidade de um componente continuar funcionando será de 65%, ou seja, a probabilidade de falha quando o componente estiver funcionando será de \[1 - p = 0.65 \Leftrightarrow p = 0.35.\]

Escolha de p e q

p <- 1 - 0.65
q <- 1 - 0.35

Algoritmo da simulação

library("tidyverse")
library("gt")

matrix_pot <- function(M, n){
  # Verificacao de pressupostos.
  # Matriz quadrada.
  if(dim(M)[1] != dim(M)[2]){
    print("Erro! A matriz nao quadrada.") # Avisa sobre o erro.
    return(M * 0) # Retorna uma matriz nula para chamar atenção.
  }
  # Com entradas nao negativas.
  if(any(M < 0)){
    return("Erro! A matriz possui alguma entrada negativa.") # Avisa sobre o erro.
    print(M * 0) # Retorna uma matriz nula para chamar atenção.
  }
  # A soma nas linhas deve ser igual a 1.
  if(any(abs(apply(M, 1, sum) - 1) > 1e-09)){
    print("Erro! Alguma linha na matriz nao soma 1.") # Avisa sobre o error.
    return(M * 0) # Retorna uma matriz nula para chamar atenção.
  }
  x <- diag(dim(M)[1])
  while(n){
    x <- x %*% M
    n <- n - 1
  }
  return(x)
}

simula_cadeia <- function(S, P, pi0, n){
  # Simula uma cadeia de Markov.
  #
  # Parâmetros de entrada:
  #
  ## S: espaco de estados finito de onde obteremos os estados indexados nas linhas e colunas da matriz de transição.
  ## P: matriz de transicao.
  ## pi0: vetor de probabilidades dos estados S.
  ## n: número de etapas a serem realizadas.

  
  # Gerar X0 com distribuicao pi0 e armazenar no vetor x.
  x <- numeric(n + 1)
  x[1] <- sample(S,size = 1, prob = pi0)
  
  # Gerar os estados das próximas etapas.
  for(i in 1:n){
    j <- match(x[i], S) # pega o índice do estado atual.
    p <- P[j,] # pega a função de transição para o estado da etapa presente.
    x[i+1] <- sample(S, size = 1, prob = p) # obter estado da proxima etapa.
  }
  return(x)
}

## Aplicação

# Já obtidos anteriormente.
p <- p # probabilidade de falha
q <- q # probabilidade de reparo

# Número de realizações
N <- 1000

P <- matrix(
  c(
    (1-q)**2, 2*q*(1-q), q**2,
    p*(1-q), p*q + (1-p)*(1-q), (1-p)*q,
    p**2, 2*p*(1-p), (1-p)**2
  ), byrow = TRUE, nrow = 3
)

# Fixar semente.
set.seed(123)

# Espaço de estados e  distribuição inicial.
S <- c(0,1,2)
pi0 <- c(0,0,1)

# Simular a cadeia.
cadeia <- simula_cadeia(S, P, pi0, N)[-1]

cadeia_tbl <- tibble(etapa = 1:length(cadeia), estado = cadeia) |>
  arrange(etapa)

tabela_resultados <- cadeia_tbl |>
  group_by(estado) |>
  count() |>
  mutate(freq = n/N) |>
  gt(groupname_col = FALSE) |>
  tab_header(paste("Simulação: número de visitas a cada estado para", N, "etapas")) |>
  cols_label(
    estado  = "Estado",
    n = "Número de visitas",
    freq = "Proporção de visitas"
  )

Resultados

• Podemos observar que a tabela com os resultados da simulação apresenta resultados compatíveis com a distribuição estacionária limite obtida de forma numérica.

Resultados

tabela_resultados

M <- P; rownames(M) <- S; colnames(M) <- S; matrix_pot(M, 999);

          0     1      2
[1,] 0.1225 0.455 0.4225
[2,] 0.1225 0.455 0.4225
[3,] 0.1225 0.455 0.4225

Resultados

• Além disso, verifica-se que a distribuição estacionária e limite obtida computacionalmente coincide com a obtida analiticamente. Com \(p = 0.35\) e \(q = 0.65\) tem-se:

\[ \begin{cases} \pi_{0} = \frac{p^{2}}{(p+q)^{2}} = \frac{0.35^{2}}{1^{2}} = 0.1225 \\ \pi_{1} = \frac{2q}{p} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{2pq}{(p+q)^{2}} = \frac{2 \cdot 0.35 \cdot 0.65}{1^{2}} = 0.4550\\ \pi_{2} = \frac{q^{2}}{p^{2}} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{q^{2}}{(p+q)^{2}} = \frac{0.65^{2}}{1^{2}} = 0.4225 \end{cases} \]

Discussão

• De acordo com os resultados da simulação anterior, a frequência de visitas aos estados em que o sistema está em operação (1 e 2) é bem maior que se comparada à frequência do estado em que o sistema está inoperante.

• Isto se deve ao fato de que especificamos uma alta probabilidade do sistema se manter em funcionamento (65%) e uma baixa probabilidade do sistema permanecer inoperante (35%).

Discussão

Confiabilidade

• Os valores de \(p\) e \(q\) impactarão diretamente na qualidade e robustez do sistema como um todo.

• Analisando-se uma máquina isolada (vide exemplo em aula de cadeia com dois estados), tem-se que a disponibilidade \(a\) de uma máquina no n-ésimo estado, isto é, \(P(X_{n} = 1)\) é dada por \(a = \frac{q}{p+q}\)

• Ou seja, em regime estacionário, o número de máquinas funcionando segue uma distribuição \(Binomial(2, a)\), e portanto o número esperado de máquinas em funcionamento é \(2a\).

Impactos de uma probabilidade reparo nula

• Em um cenário hipotético em que \(q = 0\) a matriz de transição da cadeia se torna:

\[ P = \begin{pmatrix} 1 & 0 & 0\\ p & (1-p) & 0\\ p² & 2p(1-p) & (1-p)² \end{pmatrix} \]

• Neste caso a cadeia é redutível, com o estado 0 absorvente e os demais estados transientes, portanto inevitavelmente as duas máquinas permanecerão quebradas em um n-ésimo tempo;

Impactos de uma probabilidade reparo nula

• Reordenando a matriz na forma canônica:

\[ P = \begin{pmatrix} 1-p & 0 & p\\ 2p(1-p) & (1-p)^{2} & p^{2}\\ 0 & 0 & 1 \end{pmatrix} \]

• Logo:

\[ Q = \begin{pmatrix} 1-p & 0 \\ 2p(1-p) & (1-p)^{2} \end{pmatrix} \]

e

\[ R = \begin{pmatrix} p\\ p^{2} \end{pmatrix} \]

Impactos de uma probabilidade reparo nula

• Calculando-se a matriz fundamental \(N = (I - Q)^{-1}\)

\[ (I-Q) = \begin{pmatrix} p & 0 \\ -2p(1-p) & p(2-p) \end{pmatrix} \]

\[ N = \begin{pmatrix} \frac{1}{p} & 0 \\ \frac{2(1-p)}{p(2-p)} & \frac{1}{p(2-p)} \end{pmatrix} \]

• Portanto o tempo médio de absorção partindo de cada esteado é dado por:

\[ R = \begin{pmatrix} \frac{1}{p}\\ \frac{3-2p}{p(2-p)} \end{pmatrix} \]

Impactos de uma probabilidade reparo nula

• Observa-se que caso o sistema inicie com apenas uma máquina funcionando, o tempo médio de falha é geométrico \(\frac{1}{p}\);

• Em ambos os casos, se \(p\) se aproxima de 1, o sistema tende a ficar inoperante em um passo;

• Podemos analisar um grafico com os tempos médios de falha total para diferentes valores de \(p\), considerando-se o estado inicial 2.

Impactos de uma probabilidade reparo nula